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ABSTRACT 


The effect of constant acceleration, prior to the establishment of a steady uniform 
flow, on some of the characteristics of the resulting time-dependent flow about a circular 
cylinder has been investigated numerically. It is shown that the occurrence of a local 
maximum drag is dependent on the parameters characterizing the non-impulsive nature 
of the ambient flow. However, the onset of the wake asymmetry and the evolution of the 
lift and drag forces during the transient period depend on the characteristics of the 
numerical perturbation used to initiate the asymmetric vortex shedding. It is concluded 
that the numerical methods can predict, to varying degrees of accuracy, the behavior of the 


symmetric state and the quasi-steady-state, but not of the intermediate state. 
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NOMENCLATURE 


A, = Acceleration parameter = 1/(S/R), 


Aj, = DNd"U/dt*)/V™" 


Cp = Drag coefficient = 2F/pDV’ 

Com = Maximum drag coefficient 

C, = 4F/(xpD’dU /dt) 

C. = Lift coefficient = 2L/pDV’ 

D = Diameter of circular cylinder = 2R 
F = Drag of in-line force per unit length 
L = Lift or traverse force per unit length 
R = Radius of the circular cylinder 

Re = Reynolds number = VD/v 

Rut = Radius of the outer boundary 

S = Displacement of fluid 


S/R = = Relative displacement of fluid 

(S/R),, = Relative fluid dispiacenient at which Cp,, occurs 

(S/R), = Relative displacement of fluid during the acceleration period 
=0.5(dU/dt)2/R = V?/(DdU/dt) 


=1/A, = 0.5(Vt,/R) 


St = Strouhal number = D/VT 

T = Period of vortex shedding 

t = Time 

t, = Time at the end of the acceleration period 


Vili 





U 


V 


= Time-dependent velocity 

= Constant velocity at the end of the acceleration period 

= Time step 

= Disturbance oscillation applied to ambient flow, (in Degrees) 
= Streamfunction perturbation applied to cylinder 

= Computational grid spacing 


= Plot grid spacing 
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I. INTRODUCTION 


Unsteady flow past bluff bodies has attracted a great deal of 
attention since a number of problems of practical importance are 
unsteady. Among the numerous theoretical, numerical, and 
experimental investigations, impulsively-started steady flow about 
a circular cylinder has occupied a prominent place partly because 
of its intrinsic interest towards the understanding of the 
evolution of separation, vortex formation, growth, and partly 
because it provided the most fundamental case for the comparison 
and validation of various numerical methods and codes. However, 
neither impulsive start nor impulsive stop is physically 
realizable. The flow must be accelerated from rest to a constant 
velocity, in a prescribed manner. This fact gives rise to a series 
of new questions such as: 

I What is the effect of the initial acceleration, prior to 
the establishment of a steady uniform flow, on the 
characteristics of the resulting time-dependent flow? 

2: Are there critical values of the governing parameters above 
or below which the flow may be regarded as aimost 
impulsively-started? 

33 How does the rate of accumulation of vorticity, as well as 
its cross-wake transfer, depend on the initial history of 
the motion? 

The purpose of this investigation is to explore some of these 
questions through the use of a numerical scheme based on finite- 


difference methods. 


II. BACKGROUND STUDIES 


There is a large volume of literature which deals with 
fluctuating forces and associated vortex shedding from bluff bodies 
subjected to steady ambient flow. The circular cylinder has 
attracted by far the greatest attention. Numerous computational 
studies have been performed on flow about circular cylinders in an 
attempt to predict the Strouhal number in steady ambient flow and 
the shape and growth of the wake region in impulsively started 
steady flow. Here, only the more recent and relatively more 
accurate examples will be cited. Ta Phuoc Loc [Ref. 1] solved the 
complete unsteady Navier-Stokes equation in vorticity/stream- 
function form using a combination of second and fourth-order 
compact finite-difference schemes. He obtained short-time 
symmetric-wake solutions at Reynolds numbers of 300, 550, and 1,000 
and achieved good agreement with flow visualization results for 
both vortex size and center position. His calculations also showed 
clearly the small secondary vortices just behind the separation 
points. 

Lecointe and Piquet [Ref. 2] used several compact schemes with 
the Navier-Stokes vorticity/stream function formulation to solve 
laminar flows around circular cylinders up to a Reynolds number of 
9500. They studied both start-up and unsteady periodic phenomena. 


The predicted wake-region shape showed good agreement with 


experimental flow visualizations. Ta Phuoc Loc and Bouard [Ref. 3] 





performed calculations at Re = 3,000 and 9,500 using a fourth-order 
finite-difference technique to solve the Poisson equation for the 
stream function and a second-order technique for the vorticity- 
transport equation. They found good agreement between their 
predictions and flow visualization. The calculations were 
confined, out of necessity, to relatively short times during which 
the wake became neither asymmetrical nor turbulent. Chamberlain 
[Ref. 4] used a second-order fast Poisson solver based on FET 
methods and found an accurate solution which agreed well with 
experiment and the previous computations. Rumsey [Ref. 5] used an 
upwind-biased implicit approximate factorization algorithm to 
calculate the impulsively-started unsteady flow over a circular 
cylinder at a Reynolds number 1200 and a Mach number of 0.3. 
Rumsey’s results were in very good agreement with the previous 
calculations and showed, predictably enough, only a= slight 
compressibility effect. 

All numerical calculations using finite-difference, finite- 
element, or vortex-element methods (see e.g., Sarpkaya and Shoaff 
[Ref. 6], van der Vegt [Ref. 7], Sarpkaya [Ref. 8], Chang and Chern 
(Ref. 9]) have assumed an impulsively-started flow. No 
computational attempt was made to investigate the effect of the 
initial acceleration, prior to the establishment of a steady 
uniform flow, on the characteristics of the resulting time- 
dependent flow. 

Several experimental investigations (Refs. 10-12] of 


impulsively-started flow around circular and rectangular cylinders 


have been carried out. Bouard and Coutanceau [Ref. 10] 
investigated the shape and growth rate of the wake region behind 
the cylinder for Reynolds number between 40 and 10,000. Sarpkaya 
[Ref. 11-12] examined the evolution of the wake region and the 
development of the lift and drag forces with time for cylinders 
between Reynolds numbers of 15,000 and 120,000. Nagata et al. 
(Ref. 13] studied the start-up flow at Reynolds number between 250 
and 1200, with the majority of the experiments performed at Re = 
1200. They gave detailed results for the time-evolution of the 
vortical region, boundary-layer parameters, and profile shapes at 
this Reynolds number. Sarpkaya and Kline [Ref. 14] examined the 
impulsively-started flow about four types of bluff bodies. 
Sarpkaya and Ihrig [Ref. 15] performed experiments and vortex- 
element analysis of impulsively~started flow about rectangular 
prisms and pointed out emphatically that other than numerical 
experiments, there is no mechanical system which is capable of 
generating a truly impulsive flow. In fact, efforts to generate 
i.apulsive or uniformly-accelerated flow at high Reynolds numbers 
may be hampered by the generation of compression and rarefaction 
waves and regions of intense cavitation (in liquids). Because of 


this reason one or more acceleration parameters such as 


du 
Aap Oe (7) 
V< 





or 








d"U 
A> 2=7)7 dt* (2) 
pn yn. 


will have to be added to the list of the parameters governing the 
phenomenon in order to account for the initial history of the fluid 
motion. The other parameters are the Reynolds number Re = VD/V and 


the relative displacement of the ambient flow, given by 


(AU) p2 
S dt vte (3) 
R vee R (Rt,) - if 
and 
C 
S20.5(v—*) +(t-t,) & for t>t, (4) 
R r R 


where U is the time-dependent velocity in the interval (0<t<t,), V 
is the constant velocity arrived at the end of the acceleration 
period, R is the radius of the cylinder, t is the time, and t, is 
the duration of the acceleration period. Introducing the 


dimensionless parameter defined by 


ety © 
teV— (5) 





and 


(6) 


and taking V = R = 1, Equations (1) through (4) may be reduced to 








g=(—)=— for t<t (7) 
2t, 45, ii 
and 
S=S,+(t-t,) for t>t, (8) 


where T,=25,. 

A systematic numerical variation of the governing parameters 
for an arbitrary U(t) is extremely difficult. Thus, to make 
progress one must begin with the simplest possible unsteadiness, 


namely, with constant dU/dt, so as to be able to incorporate 


progressively more complex variations of velocity with time. 





IIIT. NUMERICAL REPRESENTATION 


A. INTRODUCTION 

The fluid is assumed to be two-dimensional, incompressible 
and viscous. The governing equations for the solution are the 
Navier-Stokes equations with the stream function and the 
vorticity as independent variables. To achieve a higher 
density of mesh points near the cylinder surface, the 
computational domain is transformed from the physical plane 
(polar coordinates) to a rectangular plane. In the 
rectangular plane, the mesh is maintained at a uniform grid 
Spacing. It is necessary to have more mesh points closer to 
the cylinder surface because in this region the gradients of 
both the vorticity and the stream function are the largest. 

A third-order in time, second-order in space, three-level 
predictor-corrector finite-difference scheme is used to solve 
the vorticity-transport equation. A Fast Poisson Solver based 
on the High Order Difference approximation with Identity 
Expansion (HODIE) and the Fast Fourier Transformation (FFT) 
provided by the IMSL mathematics library is used to solve for 


the stream function. 


B. GOVERNING EQUATIONS IN THE PHYSICAL DOMAIN 
Here only a brief description of the computational 


method is presented. A more in depth description is given by 


Wang [Ref. 16]. The unsteady Navier-Stokes equations in the 
polar coordinates, as defined by the vorticity transport 


equation and the vorticity/stream-function equation are, 


a2 -2 (2 (wad) -¥ (wd) ] whe (9) 
and 
Ve w= (10) 
where 


o and YW are the vorticity and the stream function, v is the 
kinematic viscosity, t is the time and, r and © are polar 
coordinates directions (see Figure 1). The velocity 


components in the r and © directions are defined by 


ye OW (12) 
r ow 
and 
ech an (13) 
ow 
8 


The boundary conditions for the physical problem are: 
(1) no slip and zero normal velocity on the surface of the 


cylinder 


on r=R (14) 


and; (2) the potential flow at infinity is defined as 


eu(r-<) siné (15) 


and 

oO = 0 at r=ec 
where U is the external flow and R is the radius of the 
cylinder. 


The ambient flow is prescribed by 


S S = ——aee 
0<3c(Z)y u(t (16) 


in which S/R is the relative displacement of the fluid and 


(dU/dt) is the constant acceleration. 


C. TRANSFORMATION FROM PHYSICAL TO COMPUTATIONAL DOMAIN 
The coordinate transformations required to go from the 


physical domain to the computational domain are: 


r= Rexp(aS) and O= an 
where R is the radius of the cylinder, a is a transformation 
parameter (see Figure 2). After non-dimensionalizing the 
vorticity-stream function equations the transformed equations 


reduce to: 





3 _ 8 We), 8 (MWa)-2yG 3 
g(E) 57 ae oat on ( ag °? mere (17) 

and 
Vy=g(b)o (18) 

where 
we, o (19) 
G2 an? 

and 

g(&) =a*exp (2a) . (20) 


The non-dimensional fluid velocity components in the € and n 


directions are given by 





(21) 


and 





v= (22) 
vote) 
with U and V defined as 
y=- OF (23) 
on 
and 
- of 24 
V aE (24) 
The boundary conditions in the computational domain are: 
—_ OWw_ 
- Jt =0 25 
Y aE (25) 
at € = 0 ;and 
P=2Usinh (af) sin(an), (26) 
and 
@=0 (27) 


at E= 0; and 
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WE nm) =W(E n+), (28) 


and 


(En) OE, n+ | (29) 


The effect of the transformation upon the mesh spacing in 
the physical is seen in Figure 1. The bar over the non- 
dimensional quantities will no longer be used for the sake of 
simplicity. All further quantities should be understood to be 


non-dimensional. 


D. DESCRIPTION OF FINITE DIFFERENCE SCHEME 

The finite difference scheme is really a combination of 
two different schemes as previously described. The use of a 
central difference scheme to spatially discretize equation 


(17) leads to 


g(E, » 8049. Yager 4, ger Us gO 4,g-a Vier Mier, 9 Vi-r (P11, 
1.j° Or 28Eé 28n (30) 
+2 ( Wii SPF Mpa, Misi g COs Cig) 
Re 6&2 $n ’ 


where i,j represent the ith node in the N) direction and the 
jth node in the &€ direction. The nodal velocities are 


represented by: 
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U, =~ eed Vans (31) 


1, 28n 
and 
Wage Vi, 5-1 32 
V5 he (32) 


For further simplification define 





ee Uy 441% 5, 4217 Ui, 5-194, 5-1 Vier, Maen, 97 Vi-1, 55-1, 5 
med |g (6) 26 26n (33) 
+2 ( 552172 O st Os 5-1 | Wie, 572M, 3 Mi-2.J) ] 
Re 6&2 bn? 


(34) 





Taking the central difference approximation for vorticity 
and applying a two geep;. these level, predictor corrector, 
finite difference scheme with a third order accuracy in time, 
temporal discretization determines the vorticity in the 
computational domain. Expanding ow"*' into a Taylor Series and 
substituting into equation (34) the equation arrived at is 


n+1 n n hae 2 
sad OO has eae n , Of: jAt Of 5At 


Fe EASE oh AR tae OER 4j. (35) 
At tte. AOt=3! Toney 
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The second and third order terms of the right hand side of 


equation (35) is discretized and becomes 


n+l 


GD, j -@5. ose 5 £7548 £7 5- ae j 3 (36) 
ve eC +O(At?). 


The predictor used is given by the second order Adams- 


Bashforth method, as shown 


At 


Os. i =W 4, 7? (32445 i) ' (37) 


where n is defined as the n* time step and the superscript * 
indicates the predicted value. The corrector is given by 


wi 52 =Wi, jt Bt Shr, j*8 fi, 5 -f7'3) . (38) 


The term f*' is eliminated using equations 


(37) and (38), 
leaving the predictor-corrector as, 


. 4 


wo} (Wy 7754, jy) + SAt 





Cree se (39) 


In order to calculate the first time step, 


the Euler forward 
method is used to find @ 


The predictor is 


(40) 


and the corrector is 


14 





w9 70}, +48 high hs4) (41) 


The non-dimensional time scale is defined by the relative 


. displacement S/R, as noted previously. 


E. CALCULATION OF VORTICITY ON THE CYLINDER WALL 
Vorticity on the cylinder wall is determined using a 


Taylor Series expansion of the stream function, 





= 0 aw, AE 
2)=<n(1)+( 24) age ( ZY 
w (2) =p (1) (3e) 2 (e222 — 


o Ee 
+( er) 24 oy TOMES) | 





Equations (18) and (19), and the boundary condition (25) gives 


(Sedge), (43) 
and 


on the surface of the cylinder. Substituting equations (43) 


and (44) into equation (42) and using the no slip boundary 
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condition, equation (14) gives 


en Ae AR oe 45 
w (1) Eg) 3911) (AE“) (45) 


F. CALCULATION OF THE LIFT AND DRAG COEFFICIENTS 

The lift and drag coefficients are determined from the 
contribution of the viscous forces tangential to the flow and 
the pressure forces acting normal to the surface of the 
cylinder. The viscous forces are calculated from T,="Ww. This 


relationship gives the total drag force as 


--["p,cos (0) Rd®- [""posin (w) Rad, (46) 


and the total lift force as 


=-f"p,sin (8) RoB+[““pocos (8) Rd. (47) 


After dividing both the total drag and total lift equations by 


(0.5 p U* D) and defining 


i (Dp) 
are re 48 
(+ pUE) re 


the drag coefficient reduces to 
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Oe hase ae hae 
Cy =. Pp,cos (0) do =|. @sin (6) dd, (49) 


and the lift coefficient is given 


c.2-4 [°"B,s1n (8) do+ 2 [wes (6) (50) 


The pressure coefficient is determined from the Navier-Stokes 
equation in terms of dimensionless vorticity. Once integrated 
with respect to 0, the equation is 


ee 4 fn, do 
BD, (8) =p, (0) +f (FE) |e. (51) 


Equation (51) is substituted into equations (49) and (50) to 


determine the numerical scheme for the total lift and drag 


coefficients, 
ony 
Cy=- =f. Uf (2) |, ,d0] cos (6) + oo 
Bein (0) }d0, 
and 
one 
C= = Uf (2) | -,G0] sin (6 is 
-weos (8) }de. 
17 


The radial derivative of the vorticity on the surface of 
the cylinder, used in the lift and drag coefficient 


calculations, is determined using a discrete pointwise 





approximation, 
0, ~30 ,+405,,-M5,5 2 
i sn i+ bi : 54 
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IV. DISCUSSION OF RESULTS 


A. INTRODUCTION 

The numerical experiments were carried out through the use 
of a VAX-2000 system and the IMSL library. It became quickly 
evident that the type of calculations performed would require 
a computer of greater capacity and speed. Nevertheless it was 
deemed necessary to proceed with the existing system partly to 
delineate the limitations of the numerical scheme, partly to 
determine the limits of the disposable parameters, and partly 
to compare the results of the physical experiments with those 
of the numerical predictions during the early stages of flow. 
Typically, a numerical experiment requires about 50,000 times 
more time than a physical experiment (50,000 sec/1 sec). 

The fundamental objective of the calculations was to 
determine the ranges of flow which can be calculated 
accurately. The extensive literature that exists on the 
impulsive flow has either concentrated on the initial 
symmetric state or on the asymmetric late time vortex 
shedding. However, no systematic attempt was made to 
determine the upper limit of the early stages, the upper and 
lower limits of the intermediate indeterminable state, or the 
lower limit of the quasi-steady-state, in terms of the 


parameters characterizing the artificial disturbance imposed 
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on the flow. It has long been recognized that the symmetric 
state becomes increasingly unstable and the flow sooner or 
later bifurcates into an asymmetric state. This bifurcation 
is not an instantaneous event but takes place rather gradually 
even if the disturbance is imposed suddenly. However, the 
interesting feature of all the numerical calculations is that 
the numerical noise and truncation errors are ever present and 
continue to work on the propensity of the flow to become 
naturally asymmetrical even though the results are far from 
being natural. Had one been able to devise a sufficiently 
accurate numerical scheme and a greater-precision computer, 
one could maintain a longer symmetric state. Evidently, the 
onset of asymmetry in calculations depends on the 
characteristics of the physical disturbances. The two types 
of disturbances used can never be made identical, but they may 
be made to mimic each other. Thus, it is the hope of the 
numerical experimenter that the early stages of an impulsively 
started flow is relatively immune to truncation errors and the 
imposed, reasonable, artificial perturbation can, therefore, 
be expected to compare with the physical experiments. 
However, once the flow becomes asymmetrical the period of 
transition into a quasi-steady-state depends, to varying 
degrees of intensity, on the parameters characterizing the 
numerical disturbance. For small perturbations, the quasi- 
steady-state may eventually be arrived at smoothly without the 


lift and drag overshooting first and then reducing to their 
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terminal values. The rather unfortunate aspect of the 
numerical dilemma is that the quasi-steady-state is not just 
a function of the characteristics of the perturbation. If it 
were, one would have conducted a series of numerical 
experiments, arrived at a fairly stable state, and would have 
concluded that the flow no longer remembers how it was started 
and how it cver became asymmetrical. Even though this is the 
ultimate goal of the numerical experiments, the effects of the 
unavoidable truncation errors are ubiquitous and continue to 
influence the entire history of the computed flow. Thus, one 
may never be able to arrive at an accurate solution. 

It is in view of the realization of the foregoing facts 
that the results reported herein dealt with extensive 
sensitivity calculations to determine the effects of the type 
and intensity of the perturbations, the grid size, the time 
increment, and the effect of the outer boundary of the 
computational domain on the numerical experiments. In view of 
resources and time limitations noted earlier calculations were 
confined to standard runs and their variations. For the sake 
of brevity and for the ease of future reference a run made 


with: 


1. Reynolds number = 1000 
2. Disturbance strength = 0.5 
3. Sinusoidal disturbance 


4. Run time = 40 


21 





5. Time step = 0.02 

6. Computational grid spacing = 1/64 
7. Row = 51R 

8. Disturbance interval, S/R = 3-5 


9. (S/R), = 10 
will henceforth be referred to as the standard run. 


B. CHARACTERIZATION OF PERTURBATIONS 

The problem associated with the assignment of a 
perturbation is not the making of suitable choices among a 
limited number of equally sound characterizing parameters, but 
rather the difficulty of choosing a reasonable one from among 
an infinite set of perturbations and applying it at the right 
time interval. Faced with this’- problem, previous 
investigators used many types of artificial disturbances. In 
fact, there are as many original disturbances as there are 
original papers. In the present study, two types of 
disturbances with varying intensities were used. The first, 
devised by Wang [Ref. 16], was used to perturb the 
streamfunction on the cylinder from S/R=3 to 5. The shape of 
the disturbance was a step function. The amplitude of the 
step (AVY) was the only parameter varied while keeping 


everything fixed, including all other flow features (Re, At, 


Rout) : 
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The second type of perturbation was to change the 
direction of the ambient flow one full sinusoidal cycle in the 
said S/R range. The amplitude of the sine wave (Aa) was the 
only free parameter. It is worth noting that the second type 
of disturbance is less shock-like and gradually returns the 
perturbed quantity to its initial state. 

Figures 3a and 3b show C, versus S/R for the first type, 
step wise disturbance of the streamfunction, with an amplitude 
of 0.0025. Clearly CG, increases at first, almost impulsively, 
and then gradually to a maximum value of about 1.25, decreases 
sharply at the end of the acceleration period ((S/R),=10), 
then begins to undergo lift-induced oscillations at S/R=18. 
The mean drag continues to increase for a number of reasons, 
the least of which may be the limited extent of the outer 
boundary and the reflection there from. The lift reaches a 
large amplitude after the shedding of the third vortex at an 
approximate Strouhal number of about 0.23. One must hasten to 
add that no special accuracy is implied in the stated Strouhal 
number. It was derived from only a few cycles of oscillations. 
It must also be clarified that the evolutions of two 
integrated quantities like lift and drag are not the only 
means to judge the sensitivity of the predictions to the 
disturbance characteristics. Nevertheless, they serve to 
accentuate the said sensitivity better than other quantities 
such aS pressure distribution, vorticity distribution, 


streamlines or vorticity contours, or velocity and pressure 
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distributions within the flow field. It is because of this 
reason, that the lift and drag plots will continue to be the 
only plots referred to during the discussion of the 
disturbance characteristics. 

Figures 4a and 4b show a repeat of the above example with 
the sole exception that the disturbance strength is doubled. 
As expected, the stronger the intensity of the disturbance the 
sooner the inception of the lift and drag oscillations. 
Otherwise, C, for S/R<15 remains essentially unchanged. 

A more extensive series of numerical experiments were 
conducted with the second type of disturbance partly because 
it was more natural and more importantly because unlike the 
first type, this disturbance does not at any time violate the 
boundary conditions. 

Figures 5a-5b through 10a-10b show the drag and lift 
coefficients for (S/R),=10 varying only the amplitudes of the 
Sinusoidal disturbance. Figures 5a and 5b are interesting in 
the sense that the strength of the disturbance is specifically 
assigned to be zero, and yet the asymmetry is perceptible for 
S/R>35 (see Figure 5b). As noted earlier, this is due to the 
inherent truncation errors, and its onset could have been 
delayed even further had one used higher order or more stable 
numerical schemes. The remainder of the Figures from 6a to 
10b show what one would normally anticipate: the asymmetry 
sets in sooner, and the oscillations in the drag and lift 


manifest themselves at S/R values closer to (S/R),. The 


24 





interesting feature of the larger disturbances is that the 
onset of instability does not start at proportionally smaller 
S/R values. In fact, the differences between the lift and 
drag curves for disturbances of Aa = 1.0 and 2.5 are rather 
negligible. This may be interpreted as the optimum 
disturbance to be used to mimic the natural disturbances. 
However, this is not quite true since the three dimensionality 
of the disturbances as well as the resulting instabilities 
encountered in nature may result in significantly different 
flow characteristics during the transition period. As far as 
the calculations to be reported herein, a disturbance 


amplitude of Aa = 0.5 was chosen. 


C. GRID SIZE, TIME-STEP, AND BOUNDARY EFFECTS 

These will be discussed in conjunction with each other 
because of the simple fact that none can be independently 
varied without violating either the stability constraints or 
the reflections from the outer boundary. 

Figures lla through lld show the drag and 1ift 
coefficients, and the pressure and vorticity distributions as 
obtained from a standard run with the exception of the run 
time equal to 120 and R,,, equal to 150R. The outer boundary 
was extended to examine the long term effect of the 
disturbances on the quasi-steady-state region of the flow. 
This required the use of a large run time and, in turn, the 


need to minimize the reflections from the outer boundary. 
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This way, it was possible to carry out the calculations to 
S/R=120. Figure lla shows that the drag coefficient exhibits 
relatively large 1lift-induced fluctuations in the region of 
S/R values from approximately 40 to 60. These fluctuations 
have twice the period of those of lift, as expected. However, 
for S/R>60 the drag oscillations revert to the same period as 
the lift oscillations, and synchronize with the vortices 
shedding from one side of the cylinder. This tends to show 
that the vortices shed from the cylinder are larger in 
strength on one side than those from the other side. This is 
most likely due to the fact that the grid size at large S/R 
becomes too coarse and the gradients of the various vorticity 
terms become increasingly inaccurate. It is the realization 
of this somewhat anticipated fact that led to the decrease of 
the grid size. 

Normally, it would have been desirable to decrease the 
grid size and at the same time, maintain the outer boundary at 
the same distance as it was for the coarser. grid. 
Unfortunately, the speed and size constraints of the computer 
have precluded the need to maintain the size of the 
computational domain as large as R,,,=150R. Instead, for an 
otherwise standard run, the grid spacing was reduced to 1/128, 
the time step was reduced to 0.01 (to maintain the stability 
of the computer code), and R,,. was reduced to 80R. The 
results of this calculation are shown in Figures 12a thcough 


12d. For the purposes of comparison the drag and lift 
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coefficients obtained with the two grids (A&=1/64, At=0.02 and 
AE=1/128, At=0.01) are shown in Figures 13a and 13b. Clearly 
for S/R<10, the results are indistinguishable. For larger 
values of S/R, there are some differences which do not 
continue to increase, indicating that both calculations are 
understandably different and relatively stable at least for 
S/R<60. Also shown in these figures are the results of a run 
with a grid size of 1/128 and At=0.005. Surprisingly enough, 
the two runs with the grid size 1/128 and At=0.01 and At=0.005 
are virtually identical and shows that At is optimum and the 
results are dictated by the grid size, as expected. 

In summary, the analysis of the sensitivity of the 
calculations to the variation of the disposable parameters has 
shown that a standard run is capable of producing sufficiently 
stable results for S/R<40 with an outer boundary of R,,,=51R. 
In the section to follow the results obtained with the 
numerical model will be compared with selected experimental 
results, obtained at much higher Reynolds numbers (40,000 to 


60,000). 


D. SOME EXPERIMENTAL RESULTS 

Sarpkaya [Ref. 17] investigated experimentally the effect 
of constant acceleration prior to the establishment of a 
steady uniform flow on some of the characteristics of the 
resulting time-dependent flow about a circular cylinder. [It 


was shown that the occurrence of a local maximum drag, the 
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onset of wake asymmetry, and the evolution of the transverse 
force are dependent on the parameters characterizing the non- 
impulsive nature of the ambient flow. The experiments were 
conducted in a vertical water tunnel at Reynolds numbers 
considerably larger than those which could possibiy be 
predicted numerically without using a turbulence closure 
model. Even for the flow with a Reynolds number of 10vV0 the 
vortices become turbulent. It is for this reason that the 
majority of the previous calculations were confined to 
Reynolds numbers in the order of 109. There are a number of 
other fundamental differences between the experiments and the 
numerical model. In computations the flow is strictly laminar 
and two dimensional. The cylinders in nature have ends. The 
three dimensionality of the flow induced partly by the 
instability of the vortices and partly by the cylinder ends 
may cause differences of varying degrees even if the 
experiments were performed at Reynolds number corresponding to 
the calculations. 

Figures 14a through 14b for (S/R) ,=1, Figs 15a through 15b 
for (S/R) ,=5, and Figures 16a throvgh 16b for (S/R),=10 show 
respectively the drag coefficient with and without the use of 
a numerical perturbation. A careful perusal of the 
corresponding figures show that for the standard disturbance 
used, the numerical results do not differ from each other. As 
noted earlier, the effect of the disturbance exhibits itself 


for S/R<40. The reason for the selection of the S/R range 
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from 0O through 15 was to confine the comparison with 
experiments through a region relatively unaffected by the 
parameters characterizing the initial disturbances. Even 
though this statement is verified within the scope of the 
numerical calculations, it may not be correct for the physical 
experiments. Since, strictly speaking, there is no impulsive 
flow (numerical or experimental) a relatively large 
acceleration over a prescribed time period will have to be 
imposed onto the flow about a large enough cylinder to achieve 
a large enough terminal Reynolds number yielding accurately 
measurable lift and drag forces. Among the numerous obstacles 
to these objectives, one particularly stands out: vibrations. 
Thus, the very early stages of the flow (the first 0.1 or 0.2 
seconds) is accompanied not by a constant acceleration but by 
an acceleration, with a finite rate of change, superimposed on 
vibrations of high frequency. Thus, comparing the results 
with experiments one must bear in mind that the Reynolds 
nMumbers are not identical and the numerical and experimental 
flows are not created in the same manner. 

Figures 17 through 19 show plots of the numerical and 
experimental drag coefficients on the same graph for (S/R) ,=1, 
(S/R) ,=5, and (S/R),=10. The experimental results presented 
were obtained by Sarpkaya [Ref. 17]. It must be emphasi-ed 
that the discussion refers to th plots shown in the same 


graph not to their comparison. 
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For (S/R),=1 the experimental data exhibits a drag 
overshoot at S/R=4 and then the drag decreases gradually to 
its steady-state value, commonly accepted in the literature. 
It is important to note that the flow is still nearly 
symmetrical below S/R=7, the drag overshoot is not merely due 
to the symmetry of the vortices or the accumulation of 
vorticity in the symmetric vortices. Leaving aside a 
comparison of their magnitudes, but concentrating on the 
physics of the evolution of the drag coefficient in Figure 17, 
one realizes that as the vortices grow symmetrically and as 
the vorticity accumulates the drag increases to a maximum. 
Then, the elongation of the vortices the downstream motion of 
the center of vorticity, plus some mutual annihilation of 
oppositely-signed vorticity along the axis of symmetry lead to 
the decrease of the drag coefficient. Thus, it is clear that 
the drag overshoot has nothing to do with the onset of 
asymmetry. It depends on only when a symmetric pair of 
vortices will acquire optimum vorticity at an optimum 
distance. Since in the experiments, these depend on the 
initial conditions of imposed acceleration it is not 
surprising that the measured and predicted drag overshoots do 
not occur at the same S/R. 

Figures 18 through 19 show the effect of the imposed 
acceleration with greater clarity primarily because of the 
duration of the acceleration and the relatively smaller 


amplitude of the initial vibrations. In both figures, the 
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experimental values of the drag coefficient are somewhat 
smaller for S/R values smaller than (S/R),. This is 
explainable in terms of the nature of the acceleration and its 
time rate of change. For S/R values larger than (S/R),, the 
measured values are larger than those calculated. Without 
attaching undue significance, one may note in passing that C, 
(steady) equal to 1.2 for Reynolds number of 10,000 to 100,000 
and C, (steady) equal to 1.0 for a Reynolds numbers of 1000 
[Ref. 18). As far as the rise period of the drag is 
concerned, calculations contemplated in the near fu.22 will 
use the instantaneous imposed acceleration as input into the 
numerical code rather than a constant acceleration, averaged 


so as to arrive at the same (S/R), values. 


E. STREAMLINES AND VORTICITY CONTOURS 

To provide visual descriptions of the evolution of the 
velocity and vorticity fields a number of plots with varying 
degree of close-ups were made. Figure 20 shows for a standard 
run, but with no disturbance, the symmetric vortex pattern at 
S/R=5, Figures 2la through 2lc show at S/R=15 the streamlines, 
determined at various grid spacings (AA) about symmetrically 
evolving vortices. Figure 2lc is of some importance in the 
sense that it shows two phenomenon. The first is physical and 
related to the so called @ phenomena near the downstream 
shoulders of the cylinder. The other is purely numerical and 


shows the instability that evolves in the interpolation scheme 
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in a fine grid (here plot grid spacing is 1/160). The said 
instability has nothing to do with the flow instability but 
may eventually lead to sufficiently large asymmetries to cause 
vortex shedding as S/R>40 (see Figure 5b). 

Figures 22a and 22b show the streamlines and vorticity 
contours at S/R=40 for a standard run but without artificial 
disturbances. There is an almost imperceptibly small asymmetry 
(see 8 and 10 O’clock positions) as anticipated on the basis 
of Fig. 5b. Figures 23a and 23b show again at S/R=40 the 
vortex shedding due to the imposed sinusoidal perturbation 
with an amplitude of Aa = 0.1. A superposition of the 
streamlines (see Figure 23a) and the vorticity lines (see 
Figure 23b) shows that the center of vorticity is not at the 
apparent vortex center, a fact which has been known for a long 
time. This is entirely due to the time dependent nature of the 
wake. Figures 24a and 24b show the streamlines and vorticity 
lines at again S/R=40 for a larger sinusoidal disturbance (Aa 
= 0.5). The position of the separation points are not easily 
identifiable in the foregoing figures because of the crowding 
of the streamlines near the shoulder of the cylinder. It would 
be preferable to plot the instantaneous velocity profiles 
along radial lines to delineate the position of the excursion 
of the separation points. 

Finally, Figures 25a through 25d show the streamlines, the 
vorticity contours, the pressure distribution, and the 


vorticity distribution at S/R=40 for a standard run, but using 
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the first type of numerical perturbation with a step amplitude 
of AY = 0.005. The first two plots are similar to, but not 
directly comparable, with those shown in Figures 24a and 24b. 
The reason for this is obviously the non-unique nature of the 
flow and its strong dependence on the characteristics of the 
perturbation imposed. This is particularly true for large S/R 
values shown in these figures. The asymmetric nature of the 
pressure and vorticity is a consequence of the asymmetric 
vortex development. The lowest pressure occurring at about 80 
degrees from the stagnation point is comparable with those 


calculated by others. 
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V. CONCLUSIONS 


The investigation reported here warranted the following 
conclusions: : 


1. Even the higher order finite difference formulations of 
the governing equations can be solved for only relatively 
small Reynolds numbers. This is partly due to stability 
and computer constraints, and in part due to the 
difficulty of specifying appropriate perturbations 
forcing the flow to bifurcate into an asymmetric quasi- 
steady-state. 


2. The imposition of various states of uniform acceleration 
leads to results which have not been previously noted by 
others. Namely, drag rises to a finite value due to the 
added mass effect as soon as the acceleration is imposed, 
then remains fairly constant for a relative displacement 
less than about two and continues to rise to the end of 
the acceleration period or until the drag overshoot 
occurs, depending on whichever comes first. 


3. For almost impulsively-started flows the drag overshoot . 
occurs near S/R=4. For (S/R),>5, this overshoot is 
obsured by the effect of acceleration and continues to 
increase the drag to values larger than the drag 
overshoot at (S/R) ,=4. 


4. The early stages of the flow, i.e. S/R<15 can be 
calculated within the limits of the accuracy of the 
computational scheme. The results are essentially 
independent of the characteristics of the perturbation 
even if they were imposed at the start of the motion. 


5. The experimentally observed drag overshoot for almost- 
impulsively started flow occurs in the range 4<S/R<5, 
depending on the noise imposed on the flow at the early 
stages of the motion. 


6. There is a range of S/R values, for both impulsively and 
non-impulsively started flows, which is not amenable to 
correct numerical simulation. Because this region depends 
on the parameters characterizing the perturbations which 
are unknown and unknowable in physical experiments. 
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There is a third region of the flow in which the 
transient state evolves into a quasi-steady-state. It is 
assumed, for all intents and purposes that the flow does 
not remember how it arrived at the quasi-steady-state. It 
is tacitly assumed that the final state does not depend 
on the disturbances even though the nonlinear coupling of 
the disturbances and truncation errors may lead to 
somewhat different steady states. The existing computer 
speed and capacity does not allow one to increase the 
computational domain to arrive at a steady state at the 
specified Reynolds number. 


It is a fitting summary that where calculations can be 
made accurately (S/R<15) experiments are contaminated 
with noise, where the experiments can be relied upon 
(S/R>15) calculations cannot be carried out at the 
corresponding Reynolds numbers, and where the symmetric 
state bifurcates into an asymmetric’ state the 
computations and experiments cannot be compared because 
the computation depends on the disturbance (15<S/R<40) 
which cannot be imitated by the numerical perturbations 
even at the smaller Reynolds numbers. Thus, future 
researchers must concentrate on conducting numerical 
experiments at the early and later stages of flow, at 
Reynolds numbers smaller than 1000 and at carrying out 
painstakingly difficult force and pressure measurements 
on cylinders at the same small Reynolds numbers. It is 
only then that it will be possible to carry out a 
meaningful comparison between the numerical and physical 
experiments. 


35 


APPENDIX 


— em 


eee at ae ~a 
a ons 





t 
a 
TA ee 


AM ais 
RANG 
\) 


WN 


\ 


han pa 
tam ne ot ee a 
> an me oe -. 


\ 





Figure 1. Grid in the Physical Domain 
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Figure 22a. Streamlines at S/R=40 with No Disturbance 
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Figure 23a. Streamlines at S/R=40 with Ao=0.1 
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Figure 23b. Vorticity Contours at S/R=40 with 
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